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ABSTRACT 

The force-free (or low inertia) limit of magnetoliydrodynamics (MHD) can be applied 
to many astrophysical objects, including black holes, neutron stars, and accretion 
disks, where the electromagnetic field is so strong that the inertia and pressure of the 
plasma can be ignored. This is difficult to achieve with the standard MHD numerical 
methods because they still have to deal with plasma inertial terms even when these 
terms are much smaller than the electromagnetic terms. Under the force free approxi- 
mation, the plasma dynamics is entirely determined by the magnetic field. The plasma 
provides the currents and charge densities required by the dynamics of electromagnetic 
fields, but these currents carry no inertia. We present a high order Godunov scheme to 
study such force-free electrodynamics. We have implemented weighted essentially non- 
oscillatory (WENO) spatial interpolations in our scheme. An exact Riemann solver 
is implemented, which requires spectral decomposition into characteristic waves. We 
advance the magnetic field with the constrained transport (CT) scheme to preserve 
the divergence free condition to machine round-off error. We apply the third order 
total variation diminishing (TVD) Runge-Kutta scheme for the temporal integration. 
The mapping from face-centered variables to volume-centered variables is carefully 
considered. Extensive testing are performed to demonstrate the ability of our scheme 
to address force-free electrodynamics correctly. We finally apply the scheme to study 
relativistic magnetically dominated tearing instabilities and neutron star magneto- 
spheres. 
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1 INTRODUCTION 

Many high energy astrophysical objects such as black hole 
magnetospheres, ultrastrongiy magnetized neutron stars 
(magnetars) and probably relativistic jets in active galactic 
nuclei and gamma ray bursts involve relativistically magnet- 
ically dominated plasma. In those situations, the magnetic 
energy density conspicuously exceeds the thermal and rest 
mass energy density of particles. Magnetic fields play crucial 
roles in the dynamics of these astrophysical scenarios. They 
drive the outflows from astrophysical black holes, neutron 
stars and the surrounding accretion disks. The magnetic dis- 
sipation can also give rise to remarkable non-thermal emis- 
sions in the these high energy astrophysical phenomena. 
Force- free electrodynamics are believed to play an important 
role in those regimes. As the magnetic energy density greatly 
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exceeds the rest mass energy, conservative magnetohydrody- 
namic simulation often crashes in such circumstances. How- 
ever, force-free electrodynamics can behave well in such ex- 
treme magnetically dominated situations for the less impor- 
tant terms, such as the inertia and pressure, are entirely ig- 
nored in the force-free electrodynamics formulation. In force 
free electrodynamics the lorentz force p e E + J x B disap- 
pears, where p e and J are charge and current densities. This 
approximation allows us to understand the field structure of 
magnetospheres without solving for the plasma dynamics. 
However it is still generally difficult to solve FFE to find 
even stationary solutions. The first solutions of this kind 
came out for the axisymmetric aligned rotator (Contopou- 
los, Kazanas & Fendt 1999, hereafter CKF). The force-free 
constraint can be cast into an elliptic pulsar equation, which 
specifies the equilibrium configuration of the magnetic flux 
as a function of the poloidal current. CKF assumed that the 
last closed field line extends to the light cylinder. This state- 
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ment was recently relaxed (Goodwin et al. 2004; Contopou- 
los 2005; Timokhin 2006). These new steady-state solutions 
differ in spin-down energy loss, structure of current sheets 
and Y-point demarcating the magnetosphere and have been 
important for our understanding of pulsar physics. But an- 
alytic solutions can not address the stability and variabil- 
ity problems due to the assumption of stationarity. McKin- 
ney (2006) performed force-free electrodynamic simulations 
and found a unique stationary solution for the axisymmet- 
ric rotating pulsar magnetospheres. This solution is similar 
to the solution found by Contopoulos et al. (1999). But the 
force-free electrodynamic simulations of Komissarov (2006) 
failed to produce such results due to the high numerical dis- 
sipation in his code. Spitkovsky (2006) and Kalapotharakos 
& Contouplous (2009) solved the force-free electrodynamic 
equations via the Finite Difference Time Domain (FDTD) 
method. Due to the wide applications of force- free electrody- 
namics, such as magnetospheres around black holes, pulsars 
and accretion disks (Goldreich & Julian 1969; Blandford & 
Znajek 1977; Contopoulos, Kazanas & Fendt 1999; Lynden- 
Bell 2003), force-free jets stabilities (Narayan et al. 2009), we 
are motivated to develop a force-free electrodynamics code, 
taking advantage of recent progress in computational fluid 
dynamics (Toro 1997). 

Recent years have seen great progress of computational 
fluid dynamics, especially, the high order Riemann-solver 
based Godunov schemes. They have many desirable fea- 
tures and are currently believed to be superior to many 
other numerical schemes for hyperbolic systems. Godunov- 
type schemes evolve the cell-centered physical variables by 
incorporating the interactions between neighboring cells. In 
the simplest case, the first order Godunov scheme (Godunov 
1959), the distribution of conserved quantity inside each grid 
cell is assumed to be a constant. The one-dimensional inter- 
action between two such distributions is the classical Rie- 
mann problem (Toro 1997). High order schemes achieve high 
order accuracy by using high order approximations to com- 
pute the interface values and temporal updates. The first 
second-order Godunov scheme, Monotonic Upwind-centered 
Scheme for Conservative Laws (MUSCL), was proposed by 
Van Leer (1979). Piecewise Parabolic Methods (PPM) are 
high order extension of MUSCL schemes, which introduce 
parabolae as the basic interpolation functions in a zone al- 
lowing for a more accurate representation of smooth spa- 
tial gradients, as well as a steeper representation of dis- 
continuities (Colella & Woodward 1984). Essentially Non- 
Oscillatory (ENO) schemes (Harten et al. 1987) provide a 
uniformly high order accurate reconstruction, which is total 
variation bounded (TVB) reconstruction and give robust so- 
lutions for flows with discontinuities. They choose only one 
smoothest stencil out of a set of possible stencils of a fixed 
length. 

Unlike Essentially Non-Oscillatory (ENO) schemes, 
Weighted Essentially Non-Oscillatory (WENO) schemes 
(Shu 1997) take a linear combination with coefficients, called 
weights, of all possible stencils of a given size. The weights in 
the linear combination add up to unity and are distributed 
in such a way that stencils that contain discontinuities get 
extremely small weight. Further, the weights are designed in 
such a way that when the function is smooth in all stencils, 
the weights become close to the optimal ones so that the re- 
sulting linear combination of the stencils gives a higher order 



approximation - the same order as the one that the larger, 
combined, stencil would give. WENO-based schemes are ad- 
vantageous because they are able to both capture shocks 
and accurately resolve complex smooth flow structure. Sev- 
eral groups have had success in developing WENO-based 
schemes in application to relativistic astrophysics (Zhang & 
MacFadyen 2006) and cosmology (Feng et al. 2004). 

The divergence-free evolution of magnetic fields is of 
vital importance to the correct force-free electrodynamics 
(Brackbill & Barnes 1980). Unfortunately, straightforward 
application of Godunov schemes can not produce good re- 
sults for electrodynamics as the divergence free constraint 
for the magnetic field can not be automatically satisfied. 
Komissarov (2004) adopted an augmented system, where 
the divergence constraint is replaced by an additional evolu- 
tionary equation designed to minimize accumulation of error 
in any one location. In this method VB may be transported 
and dissipated like other variables. His method obtained un- 
physical fast reconnection in the current sheet that led to 
closed field line far beyond the light cylinder (Komissarov 
2006). 

Alternatively, we tried a different method, called con- 
strained transport (CT), to satisfy the divergence free con- 
straint. This method is based on a staggered configuration of 
the magnetic field and the electric field. If the initial mag- 
netic field has zero divergence, then every time step will 
maintain that to the accuracy of machine round-off error. 
To construct a numerical scheme based on the CT method 
and Godunov method, it is important that the volume- 
centered variables and the face-centered variables be cou- 
pled in a consistent manner. Usually the arithmetic average 
of face-centered values is used and second order accuracy 
can be achieved (Toth 2000). Recently, Li (2008) proposed 
a new third order method to map the face-centered vari- 
ables to the volume-centered variables. In this paper, we 
would like to make use of both the high-order Godunov 
scheme and the new mapping method for the CT method 
and construct a WENO-based Godunov-type scheme with 
constrained transport for the study of force-free electrody- 
namics. 

This paper is structured as follows: In § 2, we summa- 
rize the basic force-free electrodynamic equations. We give a 
basic description of the algorithms in § 3. Various numerical 
tests are given in § 4 and applications to physical models 
are discussed in § 5. 



2 EVOLUTION EQUATIONS OF FORCE-FREE 
ELECTRODYNAMICS 

When the plasma inertia and pressure are small, the balance 
of forces on the plasma is p e E + JxB = 0, where p e and J 
are the charge and current densities. The Maxwell equations 
together with this force-free constraint read (Gruzinov 1999; 
Blandford 2002): 

— + VxE = 0, (1) 
if = VxB-J, (2) 

^^ + [B ' (VX V (VXE)1 b - (3) 
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where p e = (V ■ E). Note that we have set the light speed 
c = 1 and 4n = 1 throughout this paper. Equation ((3} is 
a prescription for the plasma current to fulfill the force-free 
constraint. The two terms in the equation have simple phys- 
ical meanings: the perpendicular current is denoted by the 
first term. The second term is the parallel component of 
the current. The current density J is essential for nonlinear 
interactions. In practical numerical simulation, the second 
parallel term of current is cumbersome to calculate because 
it needs the interpolation of both fields and field derivatives. 
As an alternative way, we update the force-free Maxwell 
equations with only the perpendicular component of the cur- 
rent, and then remove the accumulated parallel component 
E|| of the electric field. This procedures achieve the same 
purpose of simulating a perfectly conducting plasma as the 
full equation ((3| in a simpler way. 



3 THE ALGORITHM OF FORCE-FREE 
ELECTRODYNAMICS 

The two dimensional force-free electrodynamic equations 
can be cast in the following compact form, 

dP 



dF dG _ 

dt dx dy 

where the primitive variables 
source terms S are 
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The two Jacobians A = dF/dP, B = <9G/<9P are 



A = 



B 



/ o 





V o 
/ o 






V i 









-1 















1 







-1 
















-1 







\ 

-1 





o / 

1 \ 






o / 



eigenvalues and normalized left and right eigenvectors of the 
matrix A are 



(5) 



where the prime denotes transposition. Similarly for the ma- 
trix B, we have 
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The z-direction flux of the Roe's Riemann solver is con- 
structed as follows (e.g., Roe 1981; Toro 1997): 



F(P 



k=l 



where a 
P + 1 



k j are the coefficients of the projection of 
- P~ ! . onto r k ■ , 1 ■ 
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^ a k,i+±J T k,i+i,3 
fc=l 



(8) 



Here r fe i+ i ■ are the fcth right eigenvector at the position of 
(x:,i,yj). P~ , . and P + , are the WENO spatial re- 

construction of the primitive variables at the position of 
(x i+ i,yj). The superscripts — and + denote left and right 
states, respectively. The one dimensional WENO spatial re- 
construction will be addressed in detail in the following sec- 
tions. Note that our multidimensional spatial reconstruc- 
tion is carried out in a dimension by dimension fashion. 
Consequently, the y-direction flux G i ] + ^ at the position 
( x i, Vj+l) can be constructed in a similar way. We may de- 
fine the following operator H(P"), 
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which will be used in the total variation diminishing (TVD) 
Runge-Kutta temporal update described in section 13.41 



3.1 Riemann Solver 

Our Godunov scheme requires that the numerical fluxes at 
the cell interfaces be constructed via the solution to the Rie- 
mann problems. The calculation of numerical fluxes needs 
the eigen-information of the Jacobi matrices A and B. The 



3.2 High-Order WENO Spatial Reconstruction 

We would apply Weighted Essentially Non-Oscillatory 
(WENO) reconstruction in our scheme (Shu 1997). Here 
we give a brief yet self-contained description of the one 
dimensional reconstruction method. Here we use v to 
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represent one arbitrary component of primitive variables 
P = (B x ,B y ,B z ,E x ,E y ,E x ). Consider one dimensional 
grid along the x axis, for one cell A; = {x i _i,x i+ i), the 
cell center is at x t . High-order schemes are built upon re- 
constructing some set of quantities from the cell center to 
the cell interface. We use the values of Vi in several grid 
cells around the grid cell in which the reconstruction is be- 
ing performed; this set of grid cells is called the stencil. 
Then, we combine the reconstructed profiles inside each of 
the grid cells to obtain the global reconstruction v(x). Now 
we can easily obtain the cell-interface representation of v i± i 
by evaluating v(x) at the locations of the interfaces. 

The interface value of v^ r }7 and the interface value of 

v^t are reconstructed as follows, 
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one dimensional reconstructions, described in this section, 
as building blocks for a multidimensional reconstruction. 
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with combination coefficients being c r j- Note here that 
c r] - = c, — The values of c r j for the case of k = 2 are 
give in Tab 1. The actual value of v7 1 are weighted sum of 
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as follows: 
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Following the above spatial reconstruction procedures, 
we could get the left and right reconstructed primitive vari- 
able values P~ , and P + , at the cell interface x, , i . Upon 

*+2 l +2 + 2 

the reconstructed primitive variables P~ i and P + t , we 

l +2 l +2 

could build the fluxes F i+ i, via the Riemann solver, at the 
cell interface x,, i . 

l+ 2 

To perform reconstructions in more than one dimen- 
sion, we use the dimension by dimension approach. It uses 



3.3 Constrained Transport 

Three types of methods are now applied to maintain the di- 
vergence free constraint in Godunov schemes. The first is to 
use a projection to clean the magnetic field of any divergence 
after each time step (e.g. Balsara & Kim 2003). The second 
is to add an evolutionary equation for the divergence to min- 
imize the accumulation of error in the computation domain 
(Komissarov 2006). The third is the constrained transport 
(CT) method, which designs the magnetic field difference 
equations to explicitly preserve the divergence free condition 
and is the method adopted here. We use the staggered mesh 
technique, which was proposed by Balsara & Spicer (1999), 
to preserve the divergence free condition of the magnetic 
field. We define the magnetic field components on the face 
centers and all the other quantities are defined at cell cen- 
ters. The method for constructing fluxes of the face-centered 
field (the fluxes are located at grid cell corners) based upon 
the fluxes of the cell-centered field (the fluxes are located 
at grid cell faces and are returned by a Riemann solver) 
proposed by Toth (2000) is adopted in our scheme. 

In the CT method, the integral form of the induction 
equation is based on area rather than volume averages. 
Starting from the differential form of the induction equa- 
tion (1), the magnetic fields can be updated via 
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is located at the cell corners. One can check that the above 
integration can maintain V ■ B = interior to a grid cell at 
time t n+1 if it was zero at time t n . 

To calculate the fluxes for cell-centered variables, the 
Godunov methods still need the magnetic field components 
that are defined at the cell centers. Usually the cell-centered 
magnetic fields are defined to be the average of the face- 
centered values, which is sufficient for second order accuracy. 
To achieve third order accuracy, we adopt the method pro- 
posed by Li (2008) to construct the cell-centered values from 
the face-centered values. The detail of the mapping from the 
face-centered values to the cell-centered values is described 
in the Appendix. 
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3.4 Runge-Kutta TVD Temporal Integration 



In order to achieve high order accuracy in time, we choose 
to update the variables using the high order TVD Runge- 
Kutta time integration (Shu & Osher 1988) to get P n+1 
from P n , which combines the first order Euler steps and 
involves prediction and correction. For example, the third 
order accuracy in time is achieved by 



P™ + AtH ( P 
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,n+l 



-p" + -p (2) + =-Atn[ p 



(20) 



(21) 



(22) 



where rl is the operator defined in the equation © of the 
section 13.11 Our explicit scheme is subject to the Courant- 
Friedrich-Levy (CFL) condition. For two dimensional prob- 
lem, the time step is determined by 



Ax 

c 



Ay 

c 



At — Ccfl x min 
We usually choose a CFL number Ccfl between 0.5-0.; 



4 NUMERICAL TESTS 

In this section, we present some one dimensional test prob- 
lems (Komissarov 2002, 2004) to validate the implementa- 
tion of this scheme for force-free electrodynamics. 
Fast wave : The initial solution is 

Bi = 1.0, B Z = E 2 = 0.0, 



1.0 
1.0 
0.7 



if x < -0.1 
(as + 0.1) if -0.1<x<0.1 
if x > 0.1 



Fig. [T] shows the results for a fast wave propagating in the 
positive direction. Both the head and tail of the wave are 
more resolved than in Komissarov (2002). 
Non-degenerate Alfven wave : In this problem, we consider 
the following initial solution, 

B{ = B' 2 = 1.0, E' 2 = 0.0, E' 3 = 1.0, 

r 1.0 if x < -0.1 

B' 3 = I 1.0 + |(a5 + 0.1) if -0.1<x<0.1 , 
{ 1.3 if x> 0.1 

where B' and E are measured in the wave frame that is 
moving relative to the grid with speed /3 = — 0.5. 
Degenerate Alfven wave : the initial solution is 

E[ = E' 2 =E' 3 = 0.0, B[ = 0, 

B' 2 = 2 cos (j>, B' 3 = 2 sin <f>, 



where 

r 0.0 

4>=\ ^(as + 0.1) 



if x < -0.1 

if -0.1<a;<0.1 

if s> 0.1 



0.9 
0.8 
0.7 
0.6 - 




Figure 1. The propagation of a fast wave. The time is at t = 
1.0. The upper three panels show the magnetic field components. 
The lower three panels show the electric field components. 



-0.9 r 
-1 '- 
-1.1 
-1.2 
-1.3 



-0.5 r 
-0.55 



-0.75 
-0.8 - 



Figure 2. The same as Fig. 1, but for the propagation of a non- 
dcgcncrate Alfven wave. The time is at t = 2.0. 



where B' and E are measured in the wave frame that is 
moving relative to the grid with speed /3 = 0.5. Figs. and 
[3] show the simulation results of the non-degenerate and de- 
generate Alfven wave. The scheme captures both fast and 
Alfven waves well. 

Three wave problem : The initial condition is 

B = (1.0, 1.5,3.5) E = (-1.0, -0.5,0.5) if x < 0, 
B = (1.0,2.0,2.3) E = (—1.5, 1.3, —0.5) if x > 0. 

Fig. [4] shows the results of the three- wave problem. The ini- 
tial discontinuity at x = evolves into two fast discontinu- 
ities and one stationary Alfven wave. 

A problem that evolves to B 2 — E 2 — > : The initial condi- 
tion is 



6 Cong Yu 




Figure 4. Three wave problem. The time is at t = 0.75. The 
initial discontinuity at x = evolves into two fast discontinuities 
and one stationary Alfven wave. 



B = (1.0, 1.0, 1.0) E = (0.0, 0.5, —0.5) if x < 0, 
B = (1.0, -1.0, -1.0) E = (0.0, 0.5, -0.5) if x > 0.2. 

In between the range ^ x ^ 0.2, a linear transition layer 
is set up to connect left and right state. Fig[5]shows that, in 
the transition layer, B 2 — E 2 decreases with time and finally 
the condition B 2 — E 2 > is violated. 
Stationary Alfven wave : The initial condition is 

B x = By = E z — 1.0 , E y = 0.0 , E x = — B z , 

r 1.0 if x < 

B z = \ 1.0 + 0.15{1 +sin[57r(s - 0.1)]} if < x < 0.2 
I 1.3 if x > 0.2 



Figure 6. The same as Fig. 1, but for the stationary Alfven wave 
problem. The time is at t = 1.0. 



Fig [5] shows that the wave is much more resolved than in 
Komissarov (2004), indicating the effective diffusion coeffi- 
cient is quite low. 

Current sheet : The initial condition is 
E = , B x = 1.0 , B z = 0.0 , 



B z = 



Bo if x < 
-B if x > 



The current sheet problems with Bo = 2 and Bo = 0.5 are 
shown in Figs. [7] and [8] The features are well resolved as in 
Komissarov (2004). 



5 PHYSICAL MODELS 

5.1 Neutron Star Magnetosphere Simulation 

In this section, we aims to get the structure of the neu- 
tron star magnetosphere via numerical simulations. Sim- 
ilar problems have been investigated by different meth- 
ods (Komissarov 2006; McKinney 2006; Spitkovsky 2006; 
Kalapotharakos & Contopoulos 2009) . We revisit this prob- 
lem in order to see the performance of the scheme for 
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Figure 7. The same as Fig 1., but for the current sheet problem 
with B Q = 2. The time is at t = 1.0. 




Figure 8. The same as Fig 1., but for the current sheet problem 
with B = 0.5. The time is at t = 1.0. 



more complicated field structures. We solve force-free elec- 
trodynamic equations in spherical coordinates on a uniform 
r x 6 — 400 x 800 grid. The initial magnetic field is assumed 
to be an axisymmetric dipole. The poloidal magnetic field is 
specified by the poloidal magnetic flux function ^ via 

Id* gj&A 
r sin V \r 89 ' dr J 



(B r ,Bo)^^—(-' 
r sm a \ r 

The poloidal magnetic flux function is 
* = ^ sin 2 6 , 



(23) 



(24) 



where n is the magnetic dipole moment. The electric field 
on the star is set to E = -(Oxr) xB to simulate a rotating 
conducting sphere. We then follow its evolution as the neu- 
tron star rotation is switched on. The radial inner and outer 



>^ 




-10 



-15 



Figure 9. Poloidal field lines of the neutron star magnetosphcrc. 
The vertical line is the light cylinder. The thick line is the field 
line that touches the light cylinder. 



boundaries are at O.IRlc and 2Rlc, where Rlc = c/Q. 
In our simulation, Q is set to 0.1 and the light cylinder is 
at Rlc = 10. The inner boundary conditions are the same 
as those proposed by McKinney (2006). The outer radial 
boundary conditions proposed by Godon (1996) are adopted 
as our outer radial boundary. Fig.|9]shows the poloidal mag- 
netic field lines of neutron star magnetosphere. We can see 
that the solution consists of a closed field line region extend- 
ing to the light cylinder and an open field line region. This 
solution is quite similar to that found by Contopoulos et 
al. (1999). Our simulation result also shows that the outer 
boundary condition performs quite well. Due to this bound- 
ary condition, we can put the outer boundary (at 2Rlc) 
relatively close to the light cylinder. Fig. shows that the 
closed magnetic field line beyond the light cylinder in Komis- 
sarov (2006) is avoided in our simulations, which means that 
our scheme has a lower dissipation. 



5.2 Tearing Instability in Relativistic 
Magnetically Dominated Plasma 

As an illustrative application of our scheme, we study the 
tearing instability in relativistic magnetically dominated 
plasma (Lyutikov 2006). In this simulation, to account for 
the non-ideal effects of resistivity, we adopt the same Ohm's 
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law as Komissarov, Barkov & Lyutikov (2007). We carry out 
two dimensional simulations in the x — y plane (all quanti- 
ties have only dependence on x and y, but no dependence on 
z). The computational domain is [— 20,^0] x [— 2/o,2/o] 3 with 
xq = 1.0 and yo = 2.0. We impose periodic boundary condi- 
tions at x = ±a;o and zero gradient conditions at y = ±j/o- 
The initial force-free equilibrium current sheet is: 

B = B taah(y/l)e x + B sech(y/l)e z , (25) 
E = -^pe 2 . (26) 

We adopt the units such that I = 0.1 and Bo = 1. The initial 
equilibrium current sheet is normal to the y-axis and its 
symmetry plane is y = 0. In our simulations, the resistivity 
•q is set to 10 -3 . The initial equilibrium state is perturbed by 
adding the following perturbations to the initial equilibrium 
states : 

b = (0,bosin6»(7T2:/a;o),0) , e = (0,0,0) . 

To get accurate solutions, we must have the resistive sub- 
layer well resolved. Far from the resistive sub-layer the reso- 
lution constraint is not that stringent. We adopt a variable 
resolution in the y direction in the same manner as Komis- 
sarov et al. (2007). In Fig. 1101 we show the time evolution of 
the current sheet with k = 0.314. The current sheet bocomes 
thinner in the middle of the computational domain and be- 
comes thicker at the x boundaries. Two large magnetic is- 
lands develop subsequently. By the end of the simulation 
a third smaller island forms in the middle. The simulation 
shows that our scheme well captures all the detailed physi- 
cal features inherent in this problem. Our simulation results 
also confirm the results by Lyutikov (2006) that the shortest 
growth time is equal to the geometric mean of the Alfven 
and resistive time-scales, which is a well-known result in the 
solar physics (Priest & Forbes 2000). 



6 CONCLUSIONS AND FUTURE WORK 

We have implemented the WENO-based staggered Godunov 
scheme to study the force-free electrodynamics. We advance 
the magnetic field with the constrained transport (CT) 
scheme to preserve the divergence free condition to machine 
round-off error. We apply the third order TVD Runge-Kutta 
scheme for the temporal integration. The mapping from 
face-centered variables to volume-centered variables is care- 
fully considered. A good variety of testings are performed 
to demonstrate the ability of our scheme to address force- 
free electrodynamics correctly. We also apply the scheme to 
study the neutron star magnetosphere. A force- free solution 
is found for the neutron star magnetosphere with a dipole 
surface field, the unphysical fast reconnection encountered 
by Komissarov (2006) is avoided. We finally investigate the 
tearing instability in the relativistic magnetically dominated 
plasma. Both primary and secondary magnetic island are 
well resolved in our simulations. Numerical results of the 
astrophysical models show that the scheme is robust and 
accurate. 

Strong field electrodynamics is proposed by Gruzinov 
(2008) recently. It is interesting to study different prescrip- 
tions of resistivity using our current scheme, which is left for 
future work. 



This paper mainly concerns with the two dimensional 
problems. Three dimensional extension of the scheme is 
straightforward. Another potential advantage is in its ex- 
tension to the adaptive mesh refinement (AMR) grid, which 
would allow us to study relevant problems with much higher 
spatial resolution and larger spatial extent. 

When the drift velocity exceeds the speed of light, the 
force-free electrodynamic approximation breaks down. Un- 
der such circumstances, the dynamical effects of matter be- 
comes important and can no longer be ignored. To enter 
this regime, we need to use the model of relativistic MHD 
proposed by Lyutikov & Uzdensky (2003) and Lyubarsky 
(2005). 
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APPENDIX A: MAPPING FROM 
FACE-CENTERED FIELD TO 
VOLUME-CENTERED FIELD 

We assume the magnetic field at the cell faces has the fol- 
lowing parabolic form, 

B x (x t _i,y) = alix^i) + a^x^^y + ^a f yy {x i _i)y 2 , (Al) 

B y (x, Vj _i) = b f (y j _i) + b f x (y j _ i )x+ -bLiy^^x 2 , (A2) 

where the superscript / designates that the coefficients are 
at the cell faces. To match the above two equations at the 
cell's four boundaries, the reconstructed polynomials inte- 
rior to one cell must take the form 
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B x (x,y) = a + a x x + a y y+ -a xx x + a xy xy 



1 2.1 2 . 1 3 / A o\ 

B y (x,y) = fo + + by 2/ + ifrizZ 2 + bxySJ/ 

+ ^ fe !/a2/ 2 + \b xxy x 2 y + ^b yyy y A . (A4) 

Matching the interior magnetic fields with the magnetic 
fields at cell faces, we have (Li 2008) 



J' L J_ rJR 



Ax 



a fL i a fR 

Ujyy -T 



ixyy — "!/!/!/ 



a fR _ (j/i 

"*yy "yy 



Ax 



H B + H T 



bl T - hi" 



Vxy — (isx 



Ay 



b fB + b fT 

w xx 1 ^xx 



(A5) 
(A6) 
(A7) 
(A8) 
(A9) 
(A10) 
(All) 
(A12) 
(A13) 
(A14) 



b fT - b fB 

T ^XX ^xx 

Oxxy — ~~(lxxx — 7 5 

Ay 

a o = \ (a f o L + a f o R ) ~ ^a xx (Ax) 2 , 

bo = l(b f T + b^)-±b yy (Ay) 2 , 

where the superscript L, R, T, B denote the values at the 
left, right, top, and bottom faces for a particular cell. Af- 
ter some manipulations, the volume-averaged cell-centered 
magnetic field is obtained by 

B x ,i,j ■= a + -!- (a xx {Ax) 2 + a yy (Ay) 2 ) , (A15) 



B y ,ij =bo + ^ (b xx (Ax) 2 + b yy (Ay) 2 ) 



(A16) 



The coefficients at the left and right cell faces can be calcu- 
lated as follows, 



a' y = 



B x (x i _i,y j+ i) - B x (x i _i,y 3 - 1 ) 

Ay 



(A17) 



B x {x i _i,y j+1 ) - 2B x (x i _i,y j ) + B x (x i _i,y j - 1 ) 

, (A18) 



T J 



a fR - 



(Ay) 2 

B x (x i _i,y j ) - —a f V y{Ay) 2 , 

B x {x i+ i,y j+ i) - B x (x i+ i,y j - 1 ) 
Ay 



(A19) 
(A20) 



a f R = B x (x i+ i, yj ) - ^a f y R (Ay) 2 



(A22) 



The coefficients at the top and bottom cell faces can be 
obtained similarly. In our simulations the equations (A15) 
and (A16) are used to calculate the cell-centered magnetic 
fields from the face-centered magnetic fields. 

This paper has been typeset from a TJ^X/ M}rjX file prepared 
by the author. 



B x (x i+ i,y j+ i) - 2B x (x i+ i,yj) + B x (x i+ i,yj-i) 
(AyJ 2 



(A21) 



